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ANALYTICAL  MODELING  OF 
REINFORCED  CONCRETE 
IN  TENSION 


ABSTRACT  A  smeared  crack  approach  to  fracture  of  concrete  in  mode 
I  was  implemented  in  the  finite  element  program  ADINA.  Nonlinear  con¬ 
crete  elements  with  tensile  cracking  were  modified  to  include  tensile 
strain  softening.  When  an  element  at  an  integration  point  cracks,  the 
stiffness  perpendicular  to  the  crack  is  reduced  to  zero  and  the  tensile 
stress  across  it  is  set  as  a  function  of  the  crack  opening.  Equilibrium 
iterations  were  implemented  to  redistribute  stress.  Two-  and  three- 
dimensional  models  of  a  single  edge  notched  beam  in  three-point  bending 
were  analyzed  and  compared  to  experimental  results  with  good  agreement. 
The  analytical  representation  of  mixed  mode  fracture  was  also  addressed. 
The  mechanisms  of  shear  transfer  across  a  crack  were  detailed,  and  the 
rough  crack  model,  relating  shear  stress  to  crack  opening,  is  presented 
with  discussions  on  orientation  of  successive  crack  pldnes,  tensorial 
invariance,  and  snap-back  phenomena.  Problems  are  identified  with 
modeling  bond  at  the  concrete/reinforcement  interface  and  its  effect  on 
crack  patterns. 
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of  the  crack  opening.  Equilibrium  iterations  were  implemented  to  redistribute 
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PURPOSE 


The  Office  of  Naval  Research  (ONR)  through  the  Naval  Civil  Engi¬ 
neering  Laboratory  (NCEL)  has  initiated  a  project  to  develop  fracture 
mechanics  methodology  for  design  application  of  reinforced  concrete 
elements  in  tensile  and  shear  stress  states.  Since  completion  of  an 
experimental  program  to  establish  the  fracture  energy  parameters  of 
plain  concrete  reported  in  Reference  1,  efforts  have  been  directed  to 
the  analytical  formulation  and  modeling  of  tensile  behavior  of  concrete. 
The  purpose  of  this  report  is  to  present  analytical  modeling  methodology 
of  mode  I  (crack  opening),  mixed  mode  crack  propagation  (shear),  and 
concrete  reinforcing  interface  behavior  (bond). 

This  report  supports  the  project  "Fatigue  and  Fracture  of  Concrete" 
in  the  ONR  6.1  Basic  Research  Program  YR023 . 03. OIB,  Structural  Modeling. 
In  addition  to  design  and  analysis  applications,  concrete  fracture  me¬ 
chanics  methodology  will  eventually  be  incorporated  into  damage  and  con¬ 
dition  assessment  process  of  existing  in  service  reinforced  concrete 
faci 1 i ti es . 

Concrete  cracking  in  tension  is  the  major  factor  contributing  to 
the  nonlinear  behavior  of  reinforced  concrete  elements.  The  modeling  of 
cracking  in  shear-critical  members  is  most  important  since  it  will 
determine  the  ultimate  resistance  and  post-failure  behavior. 

Crack  propagation  is  facilitated  when  the  material  is  in  a  state  of 
plane  strain.  The  material  of  thick  members  is  in  a  state  close  to 
plane  strain  in  the  interior  and  in  a  state  of  plane  stress  along  the 
edges.  An  accurate  representation  of  the  cracking  will  be  obtained  only 
if  three-dimensional  effects  are  considered.  In  the  first  part  of  this 
report  a  smeared  crack  representation  is  formulated  and  implemented  in 
the  two-  and  three-dimensional  nonlinear  concrete  elements  of  the  com¬ 
puter  code  ADINA  (Ref  2).  Experimental  results  from  tests  carried  out 
on  single  edge  notched  beams  show  good  agreement  when  compared  to  the 
analytical  predictions. 

Fracture  of  concrete  may  occur  in  three  ways:  Mode  I  (opening),  II 
(shearing)  and  III  (tearing).  Although  pure  modes  may  be  encountered, 
mixed  mode  propagation  is  more  likely.  In  the  second  part  of  this  re¬ 
port  existing  models  of  mode  II  and  mixed  mode  constitutive  relations 
are  evaluated.  Mixed  mode  crack  propagation  involves  considering  the 
transfer  of  shear  forces  across  cracks.  As  a  consequence,  successive 
crack  planes  may  form;  these  need  not  be  orthogonal  to  each  other.  This 
report  also  addresses  the  analysis  of  shear  transfer  across  cracks. 

Crack  distribution  is  greatly  affected  by  the  reinforcement- 
concrete  interface  behavior.  Proper  modeling  of  the  bond  stress-slip 
relationship  is  needed  for  an  accurate  prediction  of  the  crack  patter" 
Bond-slip  is  addressed  in  the  third  part  of  the  report. 

The  modifications  implemented  in  ADINA  in  each  case  have  been  com¬ 
piled  in  Appendixes  A,  B,  and  C. 
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PART  1.  THREE-DIMENSIONAL  MODE  I  CRACK  PROPAGATION 


INTRODUCTION 

Finite  element  modeling  of  concrete  fracture  mechanics  is  a  ver¬ 
satile  tool  for  analysis.  The  implementation  of  a  nonlinear  discrete 
crack  model  has  already  been  evaluated  (Ref  1).  However,  discrete  crack 
models  present  the  difficulty  of  varying  mesh  topology  to  represent  the 
crack  advance.  Modeling  the  crack  advance  is  further  complicated  if  the 
analysis  is  three-dimensional. 

A  smeared  crack  or  crack  band  approach  was  implemented  in  the  fi¬ 
nite  element  program  ADINA  (Ref  2).  The  fracture  zone  was  modeled  as  a 
band  of  uniformly  distributed  parallel  microcracks  having  a  blunt  front. 
This  concept  was  pioneered  by  Rashid  (Ref  3),  developed  by  Bazant  et  al. 
(Ref  4  through  7),  and  is  also  known  as  the  Crack  Band  Model  (CBM). 

After  implementation  of  the  CBM  approach  into  the  two-  and  three- 
dimensional  (2-D  and  3-D)  nonlinear  concrete  elements  of  ADINA,  the 
performance  of  the  models  was  evaluated  by  analyzing  the  results  of  the 
test  series  reported  in  Reference  1. 


THE  SMEARED  CRACK  APPROACH 

Numerous  experiments  conducted  on  tensile  specimens  have  shown  that 
after  crack  formation,  tensile  stresses  are  transferred  across  the  crack 
and  their  magnitude  decreases  with  crack  opening.  Two  stress-versus- 
crack  width  relationships  have  been  formulated  that  are  suitable  for 
finite  element  applications.  The  Fictitious  Crack  Model  (FCM)  (Ref  8) 
represents  a  single  crack  development  by  separating  the  elements  via 
introduction  of  new  nodes  and  imposes  nodal  forces  equivalent  to  the 
transferred  stresses.  The  Crack  Band  Model  transforms  the  crack  open¬ 
ing,  w,  into  strain,  e,  dividing  it  by  the  element  width  and  obtaining 
an  equivalent  stress  versus  strain  (o  -  e)  relationship. 

The  CBM  approach  assumes  that  cracked  elements  show  a  strain  soft¬ 
ening  behavior;  i.e.,  the  element  stiffness  is  negative.  This  leads  to 
a  stiffness  matrix  not  definitely  positive,  which  would  make  the  solu¬ 
tion  of  finite  element  equations  difficult  and  result  in  large  errors. 
Strain  softening  in  ADINA  is  included  for  concrete  elements  subjected  to 
stresses  beyond  the  maximum  compressive  stress.  Upon  reaching  maximum 
compressive  stress,  zero  (very  small)  stiffness  is  assigned  coupled 
with  isotropic  conditions,  and  the  stress  increments  are  computed  from  a 
uniaxial  stress-versus  strain  law.  A  similar  approach  can  be  imple¬ 
mented  in  tension.  Upon  cracking  at  an  integration  point,  a  zero  stiff¬ 
ness  can  be  assigned  across  the  crack  and  the  stress  increments  can  be 
derived  from  an  equivalent,  empirical  stress-versus-strain  law.  The 
tensile  model  becomes  orthotropic  (Ref  9).  Iterations  are  then  required 
to  satisfy  equilibrium. 
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While  transforming  the  stress-versus-crack-width  relationship  (o 
-  w)  to  stress-versus-strai n  (o  -  e),  the  element  width  is  included,  re¬ 
sulting  in  a  solution  dependent  on  element  size.  In  order  to  circumvent 
this  condition,  Bazant  and  Cedolin  (Ref  5)  suggested  linking  the  rela¬ 
tionship  to  the  fracture  energy,  G,,  forcing  the  empirical  o  -  e  law  to 
veri fy 

Gf  =  h  I  0  dE  (1) 

where  h  is  the  band  width,  e  the  strain  and  e  the  strain  beyond  which 
no  stress  is  transferred.  It  is  implied  that^the  fracture  energy  (the 
energy  needed  to  create  a  unit  fracture  area  along  the  crack  path)  is 
uniformly  distributed  across  the  width  of  the  fracture  zone  (band 
width) . 

Equation  (1)  is  verified  if 

o/f^  =  f(w/w^)  (2) 

which  satisfies 


Gf  =  I  0  dw  =  I  J  o/f^  d(w/w^) 


(3) 


whp  -o 


and 


w  =  crack  width  or  crack  opening 

w^  =  crack  width  beyond  which  no  stress  is  transferred 
f ^  =  tensi 1 e  strength 


w/w 


=  E 


o/E 


w  /w 

0  c 

w/w 


=  (e-o/E)/e  =  e/e  -  (o/fJff./E  E) 

0  0  t  t  0 

=  e/e  -  (n/f  )(e  /e  ) 

0  '  t^'  P  0 


where 


=  element  width 

E  =  Young's  modulus  of  concrete 

Ep  -  the  strain  corresponding  to  tensile  strength 

Substituting  Equation  (4)  into  Equation  (2)  yields: 

n/f^  ^  ffE/E  .o/f J 
t  o  t ' 


(4) 
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In  general  it  is  difficult  to  explicitly  obtain  stress. 

In  Reference  1,  different  o  -  w  laws  were  evaluated.  A  more  recent 
one  has  been  presented  by  Cornel issen  et  al.  (Ref  10)  and  used  by  others 
(Ref  11  and  12).  This  nonlinear  relationship  and  one  describing  linear 
softening  are  described  in  following  sections  and  have  been  implemented 
in  ADINA. 


STRESS  -  STRAIN  LAWS 
Linear  Softening 

In  spite  of  numerous  tests  on  tensile  specimens  (Ref  ID,  13,  and 
14)  showing  a  highly  nonlinear  softening,  linear  softening  is  sometimes 
used  for  computational  simplicity.  Typically  the  stress  declines 
sharply  upon  crack  initiation,  up  to  an  opening  of  approximately  15  x 
10  mm  beyond  which  the  decline  is  not  pronounced.  The  importance  of 
the  type  of  relationship  used  was  pointed  out  in  Reference  1  and  will 
also  be  demonstrated  in  the  present  study. 

A  constant  negative  softening  modulus,  E^,  can  be  defined  for 

linear  softening  and  c^c  as 

P 


2G. 


w  f. 
c  t 


(5) 


and  the  o  -  e  law  are  shown  in  Figure  la. 


Nonlinear  Softening 

The  nonlinear  o  -  w  relationship  defined  in  Reference  10  is  shown 
in  a  nond i men s i ona  1  form  in  Figure  lb.  It  is  an  empirical  formula  de¬ 
rived  by  curve  fitting  the  results  of  tensile  tests. 

o/f^  ~  { l  +  (C^w/w^)^ie^  G2w/w^)  _  yy/yy^(  j+c^^)e^  ^2^ 

where  0^  =  3  and  =  6.93. 


(6) 


From  Eguation  (3)  w  can  be  found  since  G.  and  f^.  are 
experimental  results  as  (Ascribed  in  the  following  sections, 
t i on  ( 6) : 


known  from 
From  Egua- 


1  i  a/f,  d(w/w  )  =  0. 19470 

;  u  t  O'' 


To  obtain  a  stress-versus-strain  law,  several  points  (w/w  ,o/f.) 
were  chosen  from  Equation  (6)  (see  Table  1)  and  transformed  into 
(r/t  ,o/f  ).  linear  interpolation  was  then  performed  between  data 
po i n?s . 


(a)  Linear  Softening 


Table  1.  Stress  -  Crack  Width  Relationship 
(from  Ref  10) 


w/w 

0 

o/f^ 

0.00 

1.0000 

0.05 

0.7082 

0.10 

0.5108 

0.15 

0.3817 

0.20 

0.2986 

0.25 

0.2446 

0.30 

0.2080 

0.40 

0. 1596 

0.60 

0.0904 

0.80 

0.0361 

1.00 

0.0000 

For  the  linear  approximation: 
I  J  o/f^  d(w/w^)  =  0.19704 


Since  w  -we  and  h  =  2w  (two  elements  were  chosen  across  the  crack 
band),  is'^fSond  from  ^ 

G,/2  =  0.19704  w  £  f, 
f  cot 


TEST  SERIES 

Twelve  single  edge  notched  beams  with  dimensions  102  mm  by  102  mm 
(4  in.  by  4  in.)  cross  section,  838  mm  (33  in.)  length  and  788  mm  (31 
in.)  span  were  tested  in  thr^e  point  bending.  The  notch-to-depth  ratio, 
a  /W,  was  0.5.  The  maximum  aggregate  si?e,  d  ,  was  9.5  mm  (3/8  in.). 
tRo  tests  were  carried  out  in  displacement  control. 

Concrete  properties  are  indicated  in  Table  2.  The  comp'^essive 
strength  was  measured  at  28  days  on  three  105  mm  diameter  by  305  mm  (6 
in.  by  12  in.)  standard  cylinders.  The  tensile  strength  was  obtained 
from  splitting  tensile  tests  conducted  at  28  days  on  six  similar  cylin¬ 
ders  At  the  initiation  of  the  compression  tests  strain  readings  at  the 
cylinders'  mid-height  yielded  the  modulus  of  elasticity. 

The  results  which  have  been  reported  in  detail  earlier  (Ref  1)  are 
summarized  in  Table  3.  The  fracture  energy  was  obtained  as  the  area 
undpr  the  load  versus  load-point  deflection  plot  divided  by  the  cross 
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Table  2.  Concrete  Properties 


Ingredient 

Amount 

cement 

279  kg/m^ 

water 

167  kg/m^ 

9.5  mm  gravel 

1062  kg/m^ 

sand 

907  kg/m^ 

Compressive  strength  at  28  days,  =  29.0  MPa 
Tensile  strength  at  28  days,  =  3.1  MPa 
Modulus  of  elasticity,  E  =  21.7  GPa 


sectional  area  at  the  notch.  Also  indicated  are  the  mid-span  displace¬ 
ments  at  peak  load,  d  ,  and  at  the  end  of  the  test,  d  (when  the  load 
carrying  capacity  of  ?he  beam  vanishes).  The  set-up  ffsed  is  shown  in 
Figure  2.  An  average  of  all  experimental  LLPD  plots  was  used  to  compare 
with  results  from  material  modeling  using  the  finite  element  method. 


FINITE  ELEMENT  ANALYSIS 

The  average  value  of  the  fracture  energy,  G^,  was  0.0763  N/mm. 
Values  obtained  showed  some  variation  with  size  and  configuration. 

Isoparametric  elements  and  Gaussian  integration  were  used.  The 
'smallest  number  possible  of  integration  points  was  used  (2x2  for  two 
dimensions,  2x2x2  for  three  dimensions)  to  compensate  for  tie  ex¬ 
cessive  stiffness  of  these  elements  (Ref  15). 

During  crack  propagation,  the  stiffness  matrix  had  to  be  reformed 
periodically  to  reflect  the  changes  of  stiffness  in  the  cracked  elements 
(Ref  15  and  16).  Stiffness  was  reformed  at  every  equilibrium  iteration 
for  each  loading  step  (a  full  Newton-Raphson  procedure  was  employed). 
This  approach  inherently  leads  to  a  step-size  dependency;  that  is  if  the 
loading  steps  are  too  large  then  too  few  reformations  will  be  performed 
(Ref  17).  Results  appeared  to  remain  practically  constant  for  displace¬ 
ment  steps  below  0.r''25  mm. 

Although  shear  transfer  across  a  crack  can  be  modeled  in  ADINA 
through  the  use  of  a  shear  retention  factor,  6,  no  shear  should  be 
present  in  the  center  section  of  the  symmetric  specimen.  To  avoid 
energy  consumption  a  very  low  value  of  (1  (0,0001)  was  adopted.  It 
should  be  noted  that  recent  research  is  drifting  away  from  a  constant 
shear  retention  factor  towards  a  shear  softening  approach  (Ref  12). 
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Table  3.  Experimental  Results 


Specimen 

(N/ifi) 

Peak  load 
(N) 

d 

(mS) 

(mfi) 

1 

72.3 

853 

0.17 

2.8 

2 

79.7 

999 

0.18 

2.4 

3 

85.6 

945 

0. 19 

2.8 

4 

70.5 

820 

0.17 

2.2 

5 

75.7 

910 

0.15 

3.1 

6 

72.4 

921 

0. 15 

2.2 

7 

83.4 

1011 

0.16 

2.5 

8 

75.3 

950 

0. 17 

2.9 

9 

68.1 

883 

0. 13 

2.4 

10 

68.6 

950 

0. 16 

2.6 

11 

84.1 

950 

0.16 

2.7 

12 

79.8 

997 

0. 15 

2.0 

Mean 

76.3 

932 

0. 16 

2.6 

Standard 

Deviation 

6. 1 

Optimum  configuration  for  a  finite  element  is  a  square  (2-D)  or  a 
cube  (3-D).  The  element  dimensions  in  the  fracture  zone  were  then 
chosen  as  10  mm  by  10  mm  (2-D)  or  10  mm  by  10  mm  by  10  mm  (3-D).  The 
crack  band  width  was  then  20  mm  which  is  2.1  times  the  maximum  aggregate 
size  d  and  is  consistent  with  Reference  6,  which  found  an  optimum  value 
of  the^  crack  band  width  to  be  around  3d  .  Finite  element  meshes  are 
shown  in  Figure  3. 

The  2-D  mesh  used  for  the  crack  bond  model  is  shown  in  Figure  3a. 
Only  half  of  the  specimen  is  discretized  due  to  symmetry.  In  the  C 
model,  only  nonlinear  elements  are  used.  In  the  LE  model,  only  the  five 
elements  ahead  of  the  notch  are  nonlinear  concrete  elements,  all  the 
others  are  linear  elastic  (LE  model). 

The  3-D  mesh  used  is  shown  in  Figure  3c.  Only  one  quarter  of  the 
specimen  is  discretized  due  to  the  double  symmetry.  Two  cases  were 
again  considered:  LE,  with  20  concrete  and  500  elastic  elements,  and  C, 
with  all  concrete  elements. 


FINITE  ELEMENT  RESULTS 
Two-Dimensional  Model 

To  compare  with  experimental  observations  (curve  EXP  of  Figure  4), 
the  LE  model  was  first  analyzed  using  a  linear  softening  (LE-LS)  and 
Cornel i ssen ' s  relationship  (LE-CS).  Both  load-load  point  deflation 
(LLPD)  responses  are  shown  on  Figure  4.  In  addition,  results  using 
linear  elastic  elements,  linear  softening  and  the  Fictitious  Crack  Model 
from  Reference  1  are  also  shown  (FCM-LE-LS).  The  mesh  used  with  the  FCM 
is  shown  in  Figure  3b.  Appendixes  A  and  B  detail  the  implementation  of 
the  FCM  and  the  two-dimensional  CBM. 

The  C  model  was  then  analyzed  using  Cornel i ssen ' s  softening,  and 
results  are  shown  with  the  corresponding  response  from  the  others 
(C-CS).  The  magnified  deformed  shape  is  shown  in  Figure  5. 

Three-Dimensional  Model 

LLPD  plots  for  both  LE  and  C  three-dimensional  models  are  shown  in 
Figure  6  along  with  experimental  response.  Only  nonlinear  softening  was 
used  and  the  corresponding  responses  are  marked  LE-CS  and  C-CS.  Numer¬ 
ical  and  graphic  results  shown  in  Figure  7  were  obtained  at  peak  load 
with  the  C-CS  model.  The  magnified  deformed  shape  of  the  cross  section 
at  the  notch  is  shown  in  Figure  7a  (due  to  symmetry  only  half  of  it  is 
actually  shown),  while  Figure  7b  presents  the  stresses  transferred 
across  the  cracked  elements  (at  the  integration  points).  Appendix  C 
details  the  implementation  of  the  three-dimensional  CBM. 


DISCUSSION 

Two-Dimensional  Model 

Initial  stiffness.  It  appears  from  Figure  4  that  the  CBM  approach 
yields  an  initial  stiffness  which  is  about  10%  lower  than  the  experi¬ 
mental  value.  First  cracking  occurs  at  a  load  of  about  600  N.  Since 
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most  of  the  elements  are  linear  elastic,  the  initial  response  of  the 
model  is  governed  primarily  by  mesh  geometry  and  element  size.  While 
the  FCM  analysis  used  a  fine  mesh  (205  elements)  with  a  zero  notch  width 
(Figure  3b),  the  CBM  mesh  is  much  coarser  (90  elements)  and  shows  a 
notch  20  mm  wide  (Figure  3a).  A  better  match  might  be  possible  if  the 
mesh  was  refined;  however,  the  element  should  not  be  made  smaller  than 
the  aggregate  size.  Another  possible  improvement  would  be  to  choose 
only  one  element  across  the  whole  crack  band,  cutting  the  notch  size  to 
10  mm.  However  the  latter  would  result  in  loss  of  model  symmetry  and 
the  entire  beam  would  have  to  be  discretized. 

Linear  elastic  elements  with  linear  softening.  Figure  4  indicates 
that  the  FCM  and  CBM  yield  essentially  equivalent  results  (curves  LE-LS 
and  FCM-LE-LS).  Both  theories  assume  the  total  energy  needed  to  frac¬ 
ture  the  specimen  to  be  distributed  uniformly  along  the  fracture  sur¬ 
face.  A  better  match  between  them  would  be  obtained  if  the  meshes  were 
equally  refined.  The  simplification  introduced  by  considering  linear 
strain  softening  yields  responses  further  from  the  experimental  data 
than  obtained  with  nonlinear  softening. 

Nonlinear  softening.  Using  linear  elastic  elements  with  Cornelis- 
sen's  nonlinear  stress  versus  displacement  relationship  yields  a  re¬ 
sponse  closer  to  the  experimental  behavior  (curve  LE-  CS).  Further 
refinement  with  exclusive  application  of  nonlinear  concrete  elements 
results  in  only  a  slight  response  change  from  the  linear  elastic  case 
(curve  C-CS). 

Three-Dimensional  Model 

Three-dimensional  isoparametric  elements  are  stiffer  than  two- 
dimensional  ones.  As  a  consequence  the  initial  behavior  is  closer  to 
the  experimental  data  as  shown  in  Figure  6.  After  cracking,  the  load 
obtained  for  each  displacement  step  is  slightly  lower  than  in  the  2-D 
case. 

The  use  of  only  nonlinear  concrete  elements  (C-CS)  increased  the 
computing  time  threefold  but  did  not  affect  substantially  the  response. 
The  predicted  peak  load  was  in  both  cases  lower  than  the  experimental 
value  (16%  for  C-CS  versus  13%  for  LE-CS).  It  is  reasonable  to  expect 
similarity  since  the  nonlinear  behavior  concentrates  near  the  fracture 
plane  and  the  crack  tip. 

Figure  7b  shows  the  stresses  transferred  across  the  cracked 
elements.  It  is  apparent  that  more  stress  is  transferred  towards  the 
beam's  free  edges.  It  is  concluded  that  the  crack  at  the  edges  does  not 
open  as  wide  and  propagates  slower  than  at  the  center  (which  typically 
occurs  with  metals).  However,  this  deviation  is  small  and  the  crack 
front  can  be  assumed  to  be  straight. 
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Figure  6.  Load  versus  load-point  deflection  plots,  3-D 
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7  Fracture  zone  at  peak  load. 


PART  2.  SHEAR  TRANSFER 


INTRODUCTION 

To  model  shear  transfer  across  a  crack.,  an  adequate  formulation  of 
the  constitutive  relations  representing  the  transferred  stresses  is 
needed.  Shear  transfer  yields  successive  crack  planes  which  need  not  be 
orthogonal  to  each  other.  Tensorial  invariance  is  then  addressed  for 
the  case  of  orthotropic  models.  Finally,  as  it  happens  in  the  case  of 
tensile  stress  transfer,  "snap-back"  and  instability  may  occur. 


SHEAR  TRANSFER 

Cracks  in  reinforced  concrete  are  able  to  transmit  large  shear 
forces.  Traditionally  this  transfer  has  been  neglected  because  of 
complexity  and  justified  on  the  assumption  that  this  would  be  a  con¬ 
servative  simplification.  In  some  cases  this  argument  is  erroneous  (Ref 
18  and  19).  If  a  shear  slip  occurs  along  the  crack,  the  crack  will  tend 
to  dilate.  If  the  crack  di latancy  is  prevented,  forces  normal  to  the 
crack  faces  will  appear.  These  will  have  to  be  compensated  by  tensile 
forces  on  the  reinforcement  across  the  crack,  increasing  the  potential 
for  failure. 

Shear  stresses  can  be  transferred  across  a  crack  in  three  different 
ways;  (1)  by  aggregate  interlock  as  a  result  of  the  roughness  of  the 
crack  faces,  (2)  by  dowel  action  or  shear  resistance  of  the  reinforce¬ 
ment  across  the  crack,  (3)  by  the  axial  tensile  force  component  in  the 
reinforcement  oblique  to  the  plane  of  cracking. 

For  members  with  low  reinforcement  and  for  small  crack  widths, 
aggregate  interlock  is  the  main  mechanism  of  shear  transfer.  Tests 
carried  out  on  beams  without  web  reinforcement  showed  that  aggregate 
interlock  accounted  for  up  to  75%  of  the  shear  transfer  (Ref  20).  Hence 
most  attention  will  be  given  to  this  first  mechanism  of  transfer. 


EXPERIMENTAL  BACKGROUND 

Numerous  tests  have  been  conducted  to  evaluate  the  contribution  of 
each  mechanism  of  shear  transfer.  To  assess  transfer  by  aggregate  in¬ 
terlock,  shear  displacements  were  imposed  on  concrete  specimens  with  a 
single  crack.  The  crack  width  was  maintained  constant  using  a  variable 
external  constraint  (Ref  20,  21,  and  22),  or  unconstrained  and  monitored 
(Ref  23  through  26).  In  other  cases  the  external  constraining  force  was 
maintained  constant  (Ref  27  and  28).  To  eliminate  the  effect  of  dowel 
action,  the  concrete  specimens  were  unreinforced  (Ref  20)  or  had  the 
reinforcement  through  oversized  ducts  next  to  the  crack  (Ref  23  through 
26). 
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Most  test  results  are  presented  as  families  of  curves  relating 
transferred  shear  stress  to  shear  slip  where  each  curve  corresponds  to  a 
crack  width  (Figure  8).  The  shear  stress  is  a  function  of  shear  slip 
and  crack  width  (and,  indirectly,  of  the  normal  stress). 


ANALYTICAL  BACKGROUND 

In  the  early  attemps  at  modeling  shear  transfer  in  finite  element 
methods,  the  shear  stiffness  of  a  cracked  element  was  taken  as: 

G  =  PG 
c 

where  G  is  the  shear  stiffness  of  the  uncracked  element  and  P  is  called 
the  shear  stiffness  reduction  factor.  This  was  implemented  in  ADINA. 
This  model  does  not  reflect  the  decrease  in  shear  transfer  capability 
when  the  crack  width  increases.  Shear  transfer  eventually  vanishes  as 
the  crack  width  approaches  the  aggregate  size. 

To  overcome  this  difficulty,  P  has  been  linked  to  the  crack  width 
(Ref  29  through  32).  For  instance  Cedolin  and  Dei  Poli  (Ref  30)  used: 

G  =  G  (  1  -  —  1  for  0<£<e 
c  \  j  c 

G  =  0  for  e>E 
c  c 


where 


t  =  strain  normal  to  the  crack 

=  value  of  E  after  which  there  is  no  aggregate  interlock 

When  the  crack  width  is  kept  constant  and  the  shear  slip  is  in¬ 
creased,  shear  stress  increases  to  a  plateau  independent  of  slip  (Figure 
8).  Two  analytical  models  which  represent  the  nonlinear  relationships 
between  shear  stress  and  slip  are  the  Rough  Crack  Model  of  Bazant  and 
Gambarova  (Ref  19),  and  the  Two-Phase  Model  of  Walraven  and  Reinhardt 
(Ref  23  and  24).  Both  models  will  be  more  consistent  with  experimental 
behavior  since  they  include  general  anisotropic  properties. 


THE  ROUGH  CRACK  MODEL 


The  constitutive  laws  of  the  Rough  Crack  Model  were  simplified  in 
Reference  33  as: 


0  =  a 

nn 


1^2  (^;;.2j  0.^25  y^nV 

a  -  X  ( 1- 
nt  o' 


26. 


a,  +  a. I  I  3 

d  ^  ,  4 

a  1  +  a^r 
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A  ^  =0  ^  2S  n^m 


Figure  8.  Shear  stress  versus  shear  slip. 
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These  are  empirical  expressions  based  on  Paulay  and  Loeber's 
results  (Ref  26).  The  following  assumptions  were  made: 

-  is  always  compressive 

-  for  6=0  and  6^  >  0  the  crack  faces  cannot  be  in  contact 

and  trierefore  o  =0 
nn 

"  if  6=0  there  is  no  crack  and  670  cannnt  be  obtained. 

Hence  6,  ►  0  when  6  »  0. 

t  n 


19 


for  constant  6.  and  increasinq  6  ,  both  o  and  o  . 

.  t  ^  n  '  nn '  '  nt ' 

decrease 

As  a  consequence,  if  B  is  the  crack  stiffness  matrix  defined  by 


6|a. 

b..  =  —1- 

.J  «  6j 


with  i  =  nn,nt 
j  =  n,t 


B  is  never  positive  definite  which  can  cause  numerical  problems  in 
finite  element  programs. 


IMPLEMENTATION  IN  ADINA 

The  transfer  of  tensile  stresses  across  a  crack  with  a  smeared 
crack  approach  and  tension  softening  behavior  resulted  in  a  negative 
stiffness  for  the  cracked  element  and  was  implemented  using  a  residual 
load  vector  to  redistribute  the  stresses  during  equilibrium  iterations 
(Ref  1).  This  approach  could  also  be  followed  to  include  the  transfer 
of  shear  stresses  by  combining  the  Rough  Crack  Model  and  the  Crack  Band 
Model,  discussed  in  Part  1.  A  similar  approach  for  mixed-mode  crack 
propagation  appeared  in  Reference  34. 


PRINCIPAL  AXES  ROTATION 

In  finite  element  implementations,  cracking  at  a  point  occurs  when 
one  of  the  principal  stresses  reaches  the  tensile  strength.  A  failure 
or  cracking  plane  is  defined  upon  cracking.  Most  computer  programs, 
including  ADINA,  keep  this  plane  constant  and  only  allow  successive 
planes  to  form  perpendicular  to  the  first  one  and  to  each  other.  How¬ 
ever,  when  shear  stress  transfer  across  the  crack  is  allowed,  successive 
crack  planes  will  generally  not  be  perpendicular  to  each  other  (Ref  17). 

In  order  to  address  this  inconsistency,  several  recent  approaches 
have  been  proposed,  which  are  discussed  below. 

The  Rotating  Crack  Model 

It  is  assumed  in  the  rotating  crack  model  that  the  cracks  are 
formed  normal  to  the  major  principal  tensile  strain  and  rotate  with  it. 
Experiments  carried  out  by  Vecchio  and  Collins  (Ref  35)  on  square, 
reinforced  concrete  panel  sections  support  the  assumption  that  the  main 
crack  formation  is  normal  to  the  major  principal  tensile  strain. 

Cope  et  al.  (Ref  36),  first  applied  a  rotating  crack  model,  using  a 
set  of  perpendicular  axes,  which  followed  the  tensile  strain  rotation  in 
a  step-wise  fashion.  Gupta  and  Akbar  (Ref  37)  improved  the  model  by 
considering  a  single  crack  which  followed  the  tensile  strain  rotation  in 
a  continuous  fashion.  This  method  has  been  used  by  several  researchers 
(Ref  38  through  41);  however,  the  rotating  crack  model  has  been  criti¬ 
cized  (Ref  9,  42,  and  43)  for  neglecting  the  previously  formed  cracks. 


The  Multiple  Crack.  Model 


An  alternative  approach  to  the  rotating  crack  model  was  formulated 
by  De  Borst  and  Nauta  (Ref  43,  44,  and  45)  and  independently  by  Riggs 
and  Powell  (Ref  46),  following  the  original  development  by  Litton  (Ref 
47).  They  assumed  that  multiple,  nonorthogonal  cracks  can  form  at  an 
integration  point.  In  this  procedure  the  total  strain  increment  is 
first  decomposed  into  a  solid  concrete  strain  increment  and  a  crack 
strain  increment.  Then,  crack  shear  and  normal  strains  are  related  to 
the  corresponding  stresses  and  an  incremental  crack  stress-strain  matrix 
is  derived.  In  a  similar  fashion  a  solid  concrete  stress-strain  matrix 
is  formed.  Finally,  matrices  of  all  cracks  are  assembled  and  an  ex¬ 
pression  for  the  total  stress-strain  matrix  is  obtained. 

The  multiple  crack  model  leads  to  excessive  formation  of  new 
cracks,  which  led  to  the  adoption  of  a  threshold  angle  that  allows  new 
cracks  to  form  only  after  the  rotation  of  principal  stresses  reaches 
that  angle  (Ref  43).  Numerical  difficulties  also  are  encountered  be¬ 
cause  the  crack  stress-strain  matrices  are  not  positive  definite. 


ORTHOTROPIC  VERSUS  ANISOTROPIC  MODELS 

Numerous  finite  element  analyses  have  been  conducted  using  in¬ 
crementally  linear  constitutive  equations  characterized  by  an  ortho¬ 
tropic  tangential  stiffness.  In  the  stress-free  state  isotropy  is 
assumed  and  is  replaced  by  stress-induced  orthotropy  when  the  tensile 
stress  reaches  the  tensile  strength.  This  scheme  is  also  used  in  ADINA. 
In  cases  where  the  principal  stresses  rotate  during  the  loading  history, 
this  model  is  not  tensorially  invariant;  i.e.,  the  predicted  response  is 
affected  by  the  initial  choice  of  axes  (Ref  9).  Dilatancy  of  the  crack 
cannot  be  represented  when  orthotropy  is  assumed.  Orthotropy  assumes  no 
relation  between  the  shear  strains  and  normal  stresses.  Invariance  is 
maintained  if  general  stress-induced  anisotropy  is  assumed  instead  (Ref 
9).  An  empirical  anisotropic  tangential  stiffness  matrix,  such  as  the 
one  derived  from  the  Rough  Crack  Model,  would  be  more  suitable. 


SNAP-BACK  AND  INSTABILITY 

A  general  load-deflection  response  demonstrating  the  snap-back 
phenomenon  is  shown  in  Figure  9  (Ref  48).  If  load  control  were  at¬ 
tempted  to  obtain  this  response,  the  path  ABDEJ  would  be  obtained  since 
load  control  assumes  a  monotonic  increment  of  the  load.  On  the  other 
hand,  displacement  control  would  yield  a  more  complete  response  fol¬ 
lowing  the  path  ABCDEFHI.  In  either  case  the  segment  FGH  representing 
the  snap-back  phenomenon  could  not  be  obtained  since  displacement  is 
incremented  monotonical ly . 

To  overcome  the  difficulty  in  obtaining  a  complete  response,  Riks 
(Ref  49)  and  Crisfield  (Ref  41,  48,  and  50)  developed  a  procedure,  known 
as  Riks'  or  arc-length  method,  using  a  constraint  equation  fixing  the 
step  size  in  the  load/deflection  space. 

Snap-back  may  occur  in  practice  when  strain  softening  is  con¬ 
sidered.  Two  simple  examples  are  shown  in  References  41  and  51  for  a 
bar  in  tension: 


DISPLACEMENT 

Figure  9.  Snap-back  phenomenon. 


For  a  stress-strain  law  given  by 
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where  the  bar  is  composed  of  m  square  elements  in  which  one  of  the 
elements  is  strain  softening  and  the  other  (m-1)  unloading,  the  bar  will 
have  an  average  strain  increment  of 
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It  is  observed  that  for  m>n  the  average  strain  in  the  post-peak  regime 
is  smaller  than  the  peak  load  strain  e  ;  thus,  a  snap-back  is  orginated. 
Riks'  method  has  been  implemented  in  WINA  (Ref  52)  and  applied  to  con¬ 
crete  cracking  (Ref  41,  48,  49,  and  53). 
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PART  3.  BOND-SLIP 


INTRODUCTION 

In  finite  element  analysis  of  reinforced  concrete,  bond-slip  be¬ 
tween  reinforcement  and  concrete  has  been  modeled  using  interface  ele¬ 
ments.  Interface  elements  often  use  empirical,  nonlinear  bond  stress- 
slip  relationships. 


INTERFACE  ELEMENTS 

The  simplest  interface  element  is  the  bond-link  element  developed 
by  Ngo  and  Scordelis  (Ref  54).  This  is  a  dimensionless  element  which 
connects  two  nodes  with  identical  coordinates.  It  can  be  viewed  as 
consisting  of  two  orthogonal  springs  between  the  two  nodes.  De  Groot  et 
al.  (Ref  55),  generated  a  more  complex  element  by  combining  the  rein¬ 
forcement  and  adjacent  concrete  into  a  finite  bond-zone  element. 
Hoshino  (Ref  56)  and  Schafer  (Ref  57)  developed  the  dimensionless  con¬ 
tact  element,  which  gives  a  continuous  connection  between  two  elements. 
A  comparison  between  these  different  models  (Ref  58)  showed  that  best 
results  are  obtained  using  contact  elements  with  quadratic  or  higher 
order  displacement  functions. 

An  isoparametric  contact  element  has  been  developed  by  Keuser, 
Mehlhorn  et  al.  (Ref  59,  60,  and  61),  which  is  compatible  with  the  two- 
and  three-dimensional  elements  of  ADINA.  This  contact  element  has  been 
programmed  in  a  modular  structure  to  facilitate  the  input  of  user- 
supplied  bond  stress-slip  relationship  data. 


BOND  MECHANISMS 

The  mechanism  of  bond  comprises  three  main  components:  chemical 
adhesion,  friction,  and  mechanical  interlock  between  bar  ribs  and  con^ 
Crete.  Initially,  for  very  small  values  of  bond  stress  of  up  to  1  N/mm'^ 
chemical  adhesion  is  the  only  resisting  mechanism  (Ref  62).  If  the  bond 
stress  is  increased,  chemical  adhesion  is  destroyed  and  replaced  by  the 
wedging  action  of  the  ribs.  This  wedging  action  originates  secondary 
internal  radial  cracks  (Ref  63),  longitudinal  cracks,  and  crushing  in 
front  of  the  ribs.  If  inadequate  confinement  is  provided,  bond  failure 
would  occur  as  soon  as  the  cracks  spread  across  the  concrete  cover  of 
the  bar.  With  proper  confinement,  the  bond  stress  reaches  a  maximum 
near  f  /3  before  decreasing  as  the  concrete  between  ribs  fails  in  shear 
and  a  /victional  type  of  behavior  ensues  as  shown  in  (Figure  10)  (Ref 
62). 
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Figure  10.  Typical  bond-slip  relationship. 


EMPIRICAL  BOND-SLIP  REUTIONSHIPS 

In  order  to  obtain  local  bond  stress-slip  relationships  for  finite 
element  modeling,  the  force  to  pull  short  lengths  of  embedded  rein¬ 
forcement  out  of  concrete  are  measured.  For  embedment  lengths  of  one  to 
five  lug  spacings,  consistent  bond  stress-slip  relationships  similar  to 
the  one  depicted  in  Figure  10  have  been  obtained  (Ref  62,  64,  and  65). 
The  use  of  longer  embedment  lengths  leads  to  a  nonuniform  distribution 
of  bond  stress  (Ref  66),  difficulty  in  measuring  local  values  without 
disturbing  them,  and  different  responses  if  the  bar  is  pulled  or  pushed 
(Ref  67). 

The  parameters  influencing  bond-slip  behavior  are:  load  history, 
confinement,  clear  bar  spacing,  bar  size  and  configuration,  concrete 
strength,  transverse  pressure  and  loading  rate.  Experimental-based, 
local  bond  stress-slip  relationships  have  been  derived  which  are 
adequate  for  confined  and  unconfined  bars  to  be  used  in  the  finite 
element  methodology  (Ref  65). 


EFFECT  OF  BOND  ON  CRACK  PATTERN 

Bond-slip  based  on  Oe  Groot's  model  (Ref  55)  has  been  implemented 
together  with  strain  softening  in  the  study  of  fracture  of  reinforced 
concrete  (Ref  68).  The  inclusion  of  bond-slip  produces  a  more  realistic 
crack  pattern,  which  is  less  diffuse  and  successfully  represents  primary 
and  secondary  cracking  as  observed  by  Goto  (Ref  63). 
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CONCLUSIONS  ANu  RECOMMENDATIONS 


Mode  I  Fracture 

A  smeared  crack,  approach  has  been  implemented  in  ADINA  for  two  and 
three  dimensional  nonlinear  concrete  elements  in  tension.  The  experi¬ 
mental  results  from  tests  on  12  single-edge  notched  beams  were  analyzed 
and  good  agreement  was  found.  In  particular  it  was  shown  that; 

-  the  FCM  and  CBM  approaches  yielded  similar  results 

-  the  bluntness  of  the  crack  front  affected  the  model's 
behavi or 

-  3-D  elements  yielded  a  stiffer  response  than  2-D  elements 

-  the  crack  front  could  be  assumed  straight. 

In  this  application  only  mode  I  fracture  occurs,  and  only 
tensile  stress  transfer  across  the  crack  needed  to  be  modeled. 

Mixed  Mode  Fracture 

Mixed  mode  fracture  including  shear  stress  transfer  is  the  general 
crack  propagation  mechanism.  A  benchmark  problem  in  mixed  mode  fracture 
was  presented  by  Arrea  and  Ingraffea  (Figure  11)  and  studied  in  Refer¬ 
ences  68  through  72.  Initial  attempts  at  modeling  the  shear  transfer 
using  a  constant  shear  retention  factor,  3,  and  ADINA  yielded  results 
with  almost  no  softening  after  peak  load  (Figure  12)  and  a  crack  pattern 
which  contrasts  with  experimental  observations  (Figure  13).  By  con¬ 
sidering  a  mode  II  fracture  energy.  Rots  and  De  Borst  successfully  pre¬ 
dicted  an  experimentally  verified  load-deflection  response.  However  the 
model's  crack  pattern  at  ultimate  residual  load  remained  fixed  and  was 
inconsistant  with  physical  observations. 

It  is  expected  that  the  consideration  of  an  adequate  shear  transfer 
model,  such  as  the  Rough  Crack  Model  that  includes  general  anisotropy, 
will  be  more  consistent  with  experimental  observations  and  measurements. 
By  considering  the  crack  dilatancy,  the  initial  crack  next  to  the  notch 
tip  will  tend  to  open  further  and  propagate  in  the  direction  indicated 
by  experiments. 

Bond-Sl ip 

It  is  proposed  to  update  the  current  version  of  ADINA  with  an  iso¬ 
parametric  contact  element  for  bond.  The  versatility  of  this  model  will 
allow  for  an  easy  updating  of  the  bond  stress-slip  relationship.  The 
empirical  stress-strain  relationship  derived  by  Eligehausen,  Popov  and 
Bertero  appears  most  complete  and  should  be  implemented  in  the  contact 
element. 
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Figure  12.  Load  versus  crack  mouth  sliding  displacement 


Figure  13.  Crack  pattern. 
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Appendix  A 

FICTITIOUS  CRACK  MODEL 


In  the  Fictitious  Crack  Model,  crack  propagation  is  accomplished  by 
releasing  successive  nodes  and  inserting  a  residual  force  between  them 
that  is  a  function  of  the  crack  opening.  In  this  example,  nodes  7  to  21 
are  released.  Since  these  nodes  belong  to  the  axis  of  symmetry  of  the 
specimen  the  coding  is  simplified.  A  bilinear  strain  softening  is  used. 


CHANGES  1 N  1  USE  R  .  F  77 

C*I  INSERT  USER 

C*  I 

C*I  •••  TO  SET  MFLA 

C*  I 
C*  I 

DO  100  1=7,21 
IF  (H  .EO.  I)  MFLAG  =  1 
1 00  CONT I NUE 
RETURN 

C*  I 


SUPPLIED  CODING  lUSER  58 

I  USER  59 

G  lUSER  60 

I  USER  61 
I  USER  62 


lUSER  63 
I  USE  R  64 


CHANGES  I N  USERSL  .  F77 

C*I  •••  insert  USER-SUPPLIED  CO 

C*  I 
C*  I 

XWC  =  DD(2)»2.0 

IF  (XWC  .  LE.  0,0)  then 

RR ( 2  )  =  0.0 

ELSE  IF  (XWC  .LE.  0.01840)  THEN 
RR(  2  )  =  -  5  00  .  0* (  4  .  2  -  1  82 . 56*XWC  ) 

ELSE  IF  (XWC  .LE.  0.09202)  THEN 
RR(  2  )  =  -  5  00  .  0* (  1  .  05  •  1  1  . 4 1 *XUC  ) 

ELSE 

R  R ( 2 )  =  0.0 
END  IF 

WRITE  (6,*)  ' M , DD ( 2  )  ,  R R ( 2  )  *  ,  M , DD ( 2  )  ,  RR ( 2  ) 

RETURN 

C  •  I 

C  *  F  I  L  E  END 
END 


G  USERS109 
U  S  E  R  S  1  1  0 
U  S  E  R  S  1  1  1 


U  SE  R  S  1  1  2 
USERS  1  1 3 
U  S  E  R  S  1  14 
U  S  E  R  S  1  1  5 
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Appendix  B 

CRACK  BAND  MODEL,  TW-DIMENSIONAL 


To  implement  a  smeared  crack  approach,  additional  data  has  to  be 
read  by  the  program.  The  user  must  provide  an  additional  card  (2-D 
solid  elements,  material  5,  card  d)  where  the  choice  of  softening  is 
indicated,  as  well  as  band  width,  soft  element  width  (usually  same  as 
band  width),  fracture  energy  and  maximum  aggregate  size  (format 
I5,4F10.0).  The  dimension  of  the  vector  CRKSTR  is  increased  to  memorize 
the  unloading  point  from  the  virgin  curve. 

CHANGES  I N  TODMFE . F77 


^  IDUAS/  0,  0,  0,18,18,  0,10,15,15,33,33,  0,  0,26,6*0/, 

COMMON  /SOFT/  I S C 00 E , UU C C , E L UU , G G F F , DD A A 

IF  (MODEL. £0.5}  R E A 0 ( I  I N ,  1  0 0 5  )  1  S COO E , UU C C , E L UU , G G F F , 0 0  A  A 
1005  FORMAT  (I5,AF10.0) 

COMMON  /SOFT/  I S C 0 0 E , U U C C , E L UU , G G F F  ,  0 0  A A 
URITE  (  6,2239  )  I S C 00 E , UU C C , E L UU , G G F F  ,  0 0  A  A 


T0DMFE93 
TDFE  A2 

TOFE  101 
TDFE1219 
MATRT2U 

MATRT244 


2239  FORMAT(/38H  (BB) 


/38H 

/38H 

/38H 

/38H 

/38H 

/38H 


CODE  FOR  TENSILE  STRESS  TRANSFER 
1  =  L I  NEAR  SOFTENING 
2=CORNEL I SSEN • S  SOFTENING 
SOFT  BAND  UIDTH  (UUCC) 

SOFT  ELEMENT  UIDTH  (ELUU) 
FRACTURE  ENERGY  (GGFF) 

MAXIMUM  AGGREGATE  SIZE  (DDAA) 


I  5  , 


F  1  0 . 5  , 
F  1  0  .  5  , 
F  1  0 . 8  , 
F  1  0  .  5  ) 


MATRT726 


CHANGES  I N  ELT2D4.F77 


I DU= 1 8*  I  TUO 

DIMENSION  PR0P(1),UA(18,1),YZ(1),NO05(1),NODS(1),TEMPV1(1) 
DO  10  1=1,18 


ELT2D438 

I CDMOD 1 6 

I CDMOD26 
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1 


CRKSTR(6),STRESS(4),STRAIN(4),C(4,4),N0DS{1),TEMPV1(1),  CDMOD  53 


2  TEMPV2(n,YZ(1),NOD5(1),WA(1),DUMWA{18) 

DO  1  1=1,18 

47  CALL  DCRACK  (  C  ,  S I G , A M G L E , M 00 E L  ,  I T Y P 2D , N UM C R K  ,  1  , 1 , C R K S T R  ) 

CALL  DCRACK  (  C  ,  S T R E S S , A N G , M 0 D £ L  ,  1 T Y P 2 D , M U H C R K ,  1  , 2 , C R K S T R  ) 
CALL  DCRACK  ( C , S T R E S S , A M G L E . M 0 D E L , I T Y P 2D , N U M C R K , 2 , 2 , C R K S T R  ) 
CALL  DCRACK  ( C , S T R E S S , A N G L E  ,  M 00 E L  ,  I  T Y P 2D , N UM C R K  ,  1  ,  2  ,  C R K S T R  ) 
CALL  DCRACK  { C , S T R E S S , A N G L E  .  M 00 E L  ,  I T Y P 2 D , M U M C R K  ,  1  , 2 , C R K S T R  ) 
CALL  DCRACK  ( C , S T R E S S , A N G P R  I  , MOD E L , M Y P 2D , N U M C R K  ,  1  , 2  ,  C R K S T R  ) 
CALL  DCRACK  ( C , S T R E S S  ,  A N G  ,  M OD E L  ,  I  T Y P 2D  .  N U M C R K  ,  2  ,  1  ,  C R K S T R  ) 

DO  210  1=1 , 18 


COMOO  54 
CDMOD  66 
CDMOD270 
C0M0D302 
CDMOD350 
C0M0D374 
CDMOD422 
CDM00427 
CDM00590 
CDM00596 


DIMENSION  STR(4),EPS(4),CRKSTR{6),SP1<1),SP31{1),SP32(1),SP33(1),  CRAKID15 


DIMENSION  C(4,4),SIG<4),D(4,4),T(4,4),0SIG(4),CRKSTR(6) 


DCRACK  9 


C  RELEASE  APPROPRIATE  STRESSES 

C 

98  NF=NUMCRK  ♦  1 

GO  TO  (140,120,110,155,100,100,100),  NF 
100  CALL  DSOF  ( 4 , S I  GP ,  F ALS T R , E P , CRKSTR , E , VNU , S  I  GMAT  , S I GM AC  ) 
I F  ( NUMCRK  -  5)  140,120,110 

110  CALL  DSOF  ( 2 , S I G P ,  F A L S T R , E P , C R K S T R , E , VN U , S I G M A T , S I G M A C  ) 
120  SIGP{3)=SIGP(3) 

CALL  DSOF  ( 1  , S I GP , F AL ST R , EP , CRKSTR , E , VNU , S I GM AT , S I GM AC  ) 
C 

C  ROTATE  STRESSES  TO  GLOBAL  AXES 


DCRAC204 
DCRAC205 
DCRAC206 
DCR AC207 
DCRAC208 
OCRAC209 
DCRAC21 0 
DCRAC21 1 
DCR AC2 1 2 
DCRAC2 1 3 
DCRAC2 1 4 
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SUBROUTINE  DSOF  ( IJ , S I G P , F A L S T R , E P , C R K S T R , E , V N U , S I G M A T , S I GM A C  ) 
IMPLICIT  DOUBLE  PRECISION  (  A-H,0-Z  ) 

COMMON  /SOFT/  I S C 00 E , WU C C , G G F F , D D A A , E L W U 
DIMENSION  SIGP(4),EP(4),CRKSTR(6),C0RN(11,3) 

IF  ( CRKSTR ( IJ  )  .  GT  .  0  .  DO  )  GOTO  5 
SIGP(IJ)=FALSTR 
RETURN 
5  CONT I NUE 
C 

DATA  (CORN(I,1),I=1,11)/0.,.05..1,.15,.2,.25,.3,.4,.6,.8,1.0/ 
DATA  (CORN(I, 2), 1=1, 11)/!., .7082, .5108, .3817, .2986, .2446, 

1  .2080  ,  .1596,  .  0904  ,  .0361  ,  0.0/ 

J  J  =  I  J 

IF  {JJ.EQ.4)  JJ=3 
KK=J J+3 
EEPP  =  EP(  I  J  ) 

IF  ( EP (  I  J )  . GT  . CRKST R ( KK  )  )  C R K S T R ( K K  )  = E P (  I  J  ) 

IF  ( EP (  IJ  ). L T  .  CRKSTR ( KK  )  )  E E P P  =  C R K S T R ( K K  ) 

ISS=ISCODE-2 
IF  (ISS)  10,20,30 
C 

1 0  CONT I NUE 

EETT=1/(1/E-(2*GGFF)/(SIGMAT**2*WWCC)) 

SIGP<IJ)=FALSTR+EETT*(EEPP-CRKSTR(JJ)) 

IF  ( EP<  I J  ).  L T  .  CRKSTR ( KK  )  )  S I G P ( I J  )  =  E P ( I J  ) / E E P P * S  I  G P <  I  J  ) 

IF  ( S  I  GP(  I  J )  . GT  .  FALSTR  )  S  I  G P <  I J  )  =  F A L S T R 
IF  < S I GP (  I  J >  .  L T . 0  .  DO  )  S  I  GP (  I  J  )  =  0  .  DO 
SIGP(3)=O.DO 
RETURN 
C 

20  CONTINUE 

EO=GGFF/(WWCC*0.19704*SIGMAT) 

DO  21  1=1,11 

CORN(I,3)=CORN(I,1)+CORN(I,2)*CRKSTR(JJ)/EO 
IF  ( E E PP/EO  .  L T  .  COR N (  I  , 3  )  )  GO  TO  22 

2 1  CONT  I  NUE 

22  AA=(CORN(I-1,2)-CORN(I,2))/(CORM(I-I,3)-CORM(I,3)) 

BB  =  CORN(  I  -1  ,2)-AA*CORN(  I  -1  ,3) 

SIGP(IJ)=FALSTR*<AA*EEPP/EO+8B) 

IF  ( E P (  I  J  )  .  L T  .  CRKS T R ( KK  )  )  S  I  G P (  I  J  )  =  E P (  I  J  ) / E E P P * S  I  G P (  I  J  ) 

IF  (  S  I  GP  (  I  J  )  .  G  T  .  F  A  L  ST  R  )  S  I  G  P  (  I  J  )  =  F  A  L  S  T  R 

IF  C S  I  GP (  I  J  )  .  L T  .  0  .  00  )  S  I  GP(  I  J  )  =  0 .00 

SIGP(3)=0.D0 
RETURN 
C 

30  CONTINUE 
RETURN 
C 

END 
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Appendix  C 

CRACK  BAND  MODEL.  THREE-DIMENSIONAL 


An  additional  card  is  needed  (3-D  solid  elements,  material  5,  card 
d)  with  the  same  information  as  for  the  2-D  case. 


CHANGES  1 M  THREDH.F77 

Change  at 
or  after: 


1  IDUAS  /  0,  0,  0,  25,25,  0,14,21,21,47,47,38,8*0/, 

COMMON  /SOFT/  I S C 00 E , WWC C , E L UW , G G F F , 0 D A A 

IF  (MODEL. EO. 5)  R E AD ( I  I N , 1  0 0 9 )  I S COO E , WWC C , E L U W , G G F F  ,  D 0  A A 
1009  FORMAT  < I5,4F10.0) 

COMMON  /SOFT/  I S C 00 E , WWC C , E L WW  ,  GG F F , 0 0 A A 
WRITE  (6,2239) 

2239  FORMAT(/38H  (BB)  CODE  FOR  TENSILE  STRESS  TRANSFER,  15, 

1  /38H  1=LINEAR  SOFTENING 

2  /38H  2=CORNEL I SSEN • S  SOFTENING 

3  /38H  SOFT  BAND  WIDTH  (WWCC)  ,F10.5, 

4  /38H  SOFT  ELEMENT  WIDTH  (ELWW)  ,F10.5, 

5  /38H  FRACTURE  ENERGY  (GGFF)  ,F10.8, 

6  /38H  MAXIMUM  AGGREGATE  SIZE  (DDAA)  ,F10.5) 


THRED 1 00 
THDFE  46 
THDFE102 
THDF1 190 
MATWRT14 
MATWR243 

MATWR596 


CHANGES  I N  ELT3D4.F77 


IDW=25*ITWO  ELT3D444 

DIMENSION  PR0P(1),WA(25,l),XYZ(1),NO09(’),NODS(1),TEMPV1(1)  ICH00316 

00101=1,25  I CMOD326 

1  CRKSTR(6),STRESS(6),STRAIN(6),C(6,6),RLMN(3,3),NODS(1),  CMOD3D54 


C-1 


1  TEMPV1(1),TEMPV2(1),XYZ(1).N009(1),UA{1),DUMWA(25) 

DO  1  1=1,25 

47  CALL  DCRAK3  < C , S 1 G , R L M N , M 0 D E L , N U M C R K  ,  1  ,  1  , C R K S T R  ) 

CALL  0CRAK3  (  C  ,  S  T  R  E  S  S  ,  R  L  M  N  ,  M  OO  E  L  ,  N  UM  C  R  1C  ,  1  ,  2  ,  C  R  K  S  T  R  ) 

CALL  DCRAK3  ( C , S T R E S S , R L M M , M 0 D E L , M U M C R K , 2 , 2  ,  C R K S T R  ) 

CALL  0CRAK3  ( C , S T R E S S , R L H N . MOD E L , NUH C R K , 1 , 2 , C R K S T R  ) 

159  CALL  DCRAIC3  (  C  ,  S  T  R  E  S  S  ,  R  L  M  N  ,  M  OD  E  L  ,  N  UMC  R  X  ,  1  ,  2  ,  C  R  K  S  T  R  ) 

CALL  DCRAK3  ( C , S T R E S S , R L M N , MO 0 E L , N UM C R X ,  1 , 2 , C R X S T R  ) 

130  CALL  DCRAX3  ( C , S 1 G , R L M N , M OD E L , M UM C R X  ,  2  ,  1  ,  C R X S T R  ) 

00  210  1=1,25 


CMOD3D55 
CM0D3D67 
CM0D3261 
CMOD3286 
CMOD3340 
CM0D3363 
CM0034 1 4 
CMOD3420 
CM0D3561 
CH0D3567 


DIMENSION  STR(4),EPS(4),CRXSTR(6),SP1(1),SP31(1),SP32(1),SP33(1),  CRAXID15 


DIMENSION  C ( 4 , 4 ) , S I G ( 4 ) , D (4 , 4  )  , T (4 , 4  )  , DS I G( 4  )  ,  CRXSTR ( 6  )  DCRACX  9 


C  RELEASE  APPROPRIATE  STRESSES 

C 

N  F  =  I  X  +  1 

GO  TO  (140,120,110,100,155),  NF 
100  CALL  0S0F3  ( 3 , S  I  GP  ,  F AL ST R  ,  EP  ,  CRXSTR 
110  SIGP<6)=SIGP(6) 

CALL  OSOF3  ( 2 , S  I  G P  ,  F A L S T R , E P  ,  C R X S T R 
120  SIGP(5)=SIGP(5) 

SIGP(4)=SIGP(4) 

CALL  OSOF3  ( 1 , S I G P , F A L S T R , E P , CR XS T R 
C 

C  ROTATE  STRESSES  TO  GLOBAL  AXES 


DCRAX165 
DCRAX 1 66 
DCR AX  1 67 
DCRAX 1 68 

E , VNU , S I  GHAT  , S I GMAC  )  DCRAX169 

DCRAX 1 70 

E , VNU , S I GMA T  , S I GM AC  )  0CRAX171 

DCR AX  1 72 
DCRAX 1 73 

E , VNU , S I GMAT , S I GMAC )  0CRAX174 

D  CR AX  1 75 
OCRAXI 76 
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SUBROUTINE  DSOF3  ( IJ , S I G P . F A L S T R , E P , C R K S T R . E , V N U , S I GM A T , S I G M A C  ) 
IMPLICIT  DOUBLE  PRECISION  (  A-H,0-Z  ) 

COMMON  /SOFT/  I S C 00 E , W W C C , G G F F . 0 D A A , E L WU 
DIMENSION  SIGP<4),EP(4),CRICSTR<6).C0RN(11,3) 

IF  ( CRKSTR (  IJ  )  . GT  .  0  .  DO  )  GOTO  5 
SIGP(IJ)=FALSTR 
RETURN 
5  CONTINUE 
C 

DATA  (CORN(I,1),I  =  1,in/0...05,.1..15,.2,.25,.3,.4,.6,.8,1.0/ 
DATA  (CORNC  I, 2),  1=1,  m/l... 7082,. 5108,. 3817,. 2986,. 2446,. 2080, 
1  . 1 596, . 0904, . 0361 , 0 . 0/ 

J  J  =  I  J 
KK  =  J  J  +  3 
EEPP  =  EP{  I  J  ) 

IF  ( EP ( IJ  )  .  GT  .  CRKST R ( KK)  )  C R K S T R { K K  )  =  E P (  I  J  ) 

IF  ( E P {  I J  ). L T . CRKSTR ( KK  )  )  E E P P  =  C R K S T R ( K K  ) 

ISS=ISCODE-2 
IF  (ISS)  10,20,30 
C 

1 0  CONT I NUE 

EETT=1/{1/E-(2*GGFF)/(SIGMAT**2*UUCC)) 

SIGP{IJ)=FALSTR+EETT»(EEPP-CRKSTR(JJ)) 

IF  ( EP (  I J )  . L T  . CRKST R ( KK  )  )  S  I  G P (  I  J  )  =  E P (  I  J  ) / E E P P • S  I  G P (  I  J  ) 

IF  ( S I GP(  I J )  . GT  .  FALSTR  )  S  I  G P <  I  J  )  =  F A L S T R 
IF  ( S I GP (  I J  )  .  LT . 0  .  DO  )  SIGP(  I J)  =  O.DO 
IF  <  I J  •  2  )  12,11,11 

11  SIGP(6)=0.00 

12  SIGP(5)=0.D0 
SIGP(4)=0.00 
RETURN 

C 

20  CONT I NUE 

EO=GGFF/(UUCC*0.19704*SIGHAT) 

DO  23  1=1,11 

CORN(I,3)=CORN(I,1)+CORN<I,2)»CRKSTR(JJ)/EO 
IF  ( EEPP/ EO  .  L T  .  COR N (  I  , 3 )  )  GO  TO  24 

23  CONTINUE 

24  AA=(CORN(I-1,2)-CORN<I,2))/<CORN(I-1,3)-CORN(I,3)) 
BB=CORN(I-1,2)-AA»CORN(I-1,3) 

SIGP(IJ)=FALSTR*(AA*EEPP/EO+B8) 

IF  ( S I GP(  I J )  . GT  .  FALSTR  )  S  I  G P  (  I  J  )  =  F A L S T R 
IF  ( S  I  GP (  I  J  )  .  L T  .  0  .  D 0  )  S  I  GP (  I  J  )  =0  .  DO 
IF  (IJ-2)  22,21,21 

21  SIGP(6)=0.D0 

22  SIGP(5)=O.DO 
SIGP(4)=0.D0 
RETURN 

C 

30  CONT  I  NUE 
RETURN 


r.-i 
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END 
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